# -*- coding: utf-8 -*-
# =======================================
#               流型转换
# =======================================

import warnings
import sys
import numpy as np
from utils.constant import Physics

# 禁用所有警告
warnings.filterwarnings("ignore")


class FlowPatternClassfier:
    """流型判断

    Args:
        rho_L (float): 单点的液相密度, 默认值为1000kg/m^3
        rho_g (float): 单点的气相密度, 默认值为1.28kg/m^3
        mu_g (float): 单点的气相粘性, 默认值为0.0000181
        mu_L (float): 单点的液相粘性, 默认值为0.000089
        eps (float): 管单点的绝对粗糙度, 默认值为0.00003m
        sigma (float): 单点气液相间的表面张力, 默认值为0.025N
        V_g (float): 单点气相速度, 默认值为0m/s
        V_L (float): 单点的液相速度, 默认值为10m/s
        g (float): 重力加速度, 默认值为9.8m/s^2
        alpha_L (float): 单点液相截面分数, 默认值为0.8
        alpha_g (float): 单点气相截面分数, 默认值为0.2
        D (float): 管道直径, 默认值为1m
        x (ndarray): 管道节点里程, 为二元数组, 默认值为[0, 0]
        h (ndarray): 管道节点高程, 为二元数组, 默认值为[0, 10]

    Returns:
        str: "StratifiedFlow", "AnnularFlow", "BubbleFlow" 或者 "SlugFlow

    Example:
    # >>> fp = FlowPatternClassfier() # 实例化
    # >>> fp.rho_L = 1000    # 属性值复写
    # >>> fp.V_g = 0
    # >>> print(fp.pattern())
    """

    # # 物性参数(default)
    # rho_L = 1000
    # rho_g = 1.28
    # mu_g = 0.0000181
    # mu_L = 0.000089
    # eps = 0.00003
    # sigma = 0.025

    # # 运动参数
    # V_g = 0
    # V_L = 10
    # g = 9.8

    # # 几何参数
    # alpha_L = 0.8
    # alpha_g = 0.2
    # D = 1
    # X = np.array([0, 0])
    # h = np.array([0, 10])

    # todo, 需要把默认值在初始化时赋值，不应该在外部动态赋值
    def __init__(self):
        self.g = Physics.g.value

    def pattern(self) -> str:
        """输出流型

        Returns:
            str : "StratifiedFlow", "AnnularFlow", "BubbleFlow" 或者 "SlugFlow
        """
        return self.judge_pattern()

    def area(self) -> float:
        """计算管道截面积

        Returns:
            float: 管道截面积
        """
        return (self.D / 2) ** 2 * np.pi

    def theta(self) -> float:
        """计算管道与水平面的夹角

        Returns:
            float: 管道与水平面的夹角
        """
        if (self.X[0] == 0) and (self.X[1] == 0):
            return np.pi / 2
        else:
            return np.arctan((self.h[1] - self.h[0]) / (self.X[1] - self.X[0]))

    def discount_V_g(self) -> float:
        """气相折算速度

        Returns:
            float: 气相折算速度
        """
        return self.V_g * self.alpha_g

    def discount_V_l(self) -> float:
        """液相折算速度

        Returns:
            float: 液相折算速度
        """
        return self.V_L * self.alpha_L

    def re_L(self) -> float:
        """计算液相雷诺数

        Returns:
            float: 液相雷诺数
        """
        return self.rho_L * self.V_L * self.D / self.mu_L

    def re_g(self) -> float:
        """计算气相雷诺数

        Returns:
            float: 气相雷诺数
        """
        return self.rho_g * self.V_g * self.D / self.mu_g

    #!==================================================================
    #!                    分层流型向非分层流的转换
    #!==================================================================
    def cal_liquid_level(self) -> float:
        """采用二分法计算截面的液位高度

        Returns:
            float: 截面液位高度
        """

        def f(hl):
            f = (
                np.arccos(1 - 2 * hl / self.D)
                - 1 / 2 * np.sin(2 * np.arccos(1 - 2 * hl / self.D))
                - np.pi * self.alpha_L
            )
            return f

        error = 1
        hl1 = 0  # f(0) <= 0
        hl2 = self.D  # f(D) >= 0
        half = 0
        if f(hl1) * f(hl2) <= 0:
            while error > 1e-6:
                half = hl2 - (hl2 - hl1) / 2  # 中间点
                val = f(half)
                if val < 0:
                    hl1 = hl1 + (hl2 - hl1) / 2
                elif val > 0:
                    hl2 = hl2 - (hl2 - hl1) / 2
                else:
                    return half
                error = np.abs(val)
            return half
        else:
            sys.exit("无法求得液位高度稳态解, 无法判断流型")

    def dimless_F(self) -> float:
        """计算无因次数F

        Returns:
            float: 计算无因次数F
        """
        F = (
            np.sqrt(self.rho_g / (self.rho_L - self.rho_g))
            * (self.discount_V_g() + 1e-6)
            / np.sqrt(self.D * self.g * np.cos(self.theta()))
        )
        return F

    def judge_layered_flow(self) -> str:
        """判断分层流向非分层流转换

        Returns:
            str: "non-StratifiedFlow"或者"StratifiedFlow"型流体
        """
        h_L_tilde = self.cal_liquid_level() / self.D
        V_g_tilde = self.V_g / (self.discount_V_g() + 1e-6)

        A_g_tilde = self.area() * self.alpha_g / self.D ** 2
        dAdH_L = (1 - (2 * h_L_tilde - 1) ** 2) ** 0.5

        dimless = (
            self.dimless_F() ** 2
            * (1 / (1 - h_L_tilde) ** 2 * V_g_tilde ** 2 * dAdH_L / A_g_tilde)
            - 1
        )

        if dimless >= 0:
            return "non-StratifiedFlow"
        else:
            return "StratifiedFlow"

    #!==================================================================
    #!                       环状流向段塞流的转换
    #!==================================================================
    def dimless_coefficients(self) -> tuple:
        """根据雷诺数确定计算无因次数的系数

        Returns:
            tuple: C_g, C_L, m, n
        """
        if self.re_g() < 2000:
            C_g = 64
            m = 1
        else:
            C_g = 0.184
            m = 0.2

        if self.re_L() < 2000:
            C_L = 64
            n = 1
        else:
            C_L = 0.184
            n = 0.2
        return C_g, C_L, m, n

    def Lockhart_Martinelli2(self) -> float:
        """计算无量纲数X^2, Lockhart-Martinelli参数的平方

        Returns:
            float: 无量纲数X^2
        """

        C_g, C_L, m, n = self.dimless_coefficients()
        numerator = (
            4
            * C_L
            / self.D
            * (self.rho_L * self.discount_V_l() * self.D / self.mu_L) ** (-n)
            * self.rho_L
            * self.discount_V_l() ** 2
            / 2
        )
        denominator = (
            4
            * C_g
            / self.D
            * (self.rho_g * self.discount_V_g() * self.D / self.mu_g) ** (-m)
            * self.rho_g
            * self.discount_V_g() ** 2
            / 2
        )

        if denominator == 0:
            return numerator / 1e-6
        else:
            return numerator / denominator

    def dimless_gravity(self) -> float:
        """计算无量纲重力参数Y

        Returns:
            float: 无量纲重力参数Y
        """
        C_g, _, m, _ = self.dimless_coefficients()
        numerator = (self.rho_L - self.rho_g) * self.g * np.sin(self.theta())
        denominator = (
            4
            * C_g
            / self.D
            * (self.rho_g * self.discount_V_g() * self.D / self.mu_g) ** (-m)
            * self.rho_g
            * self.discount_V_g() ** 2
            / 2
        )
        if denominator == 0:
            return numerator / 1e-6
        else:
            return numerator / denominator

    def steady_alpha_L(self) -> float:
        """稳态求解平衡态持液率, 前提是持液率在[0, 1]区间有且仅有唯一解

        Returns:
            float: 平衡态持液率, 有解则返回持液率, 无解则返回错误码-1
        """
        error = 1
        X2 = self.Lockhart_Martinelli2()
        Y = self.dimless_gravity()

        def f(alpha_L):
            f = (
                (1 + 75 * alpha_L) * alpha_L**2 / (1 - alpha_L) ** 2.5
                - X2
                - Y * alpha_L**3
            )
            return f

        hl1 = 1e-6  # f(0) <= 0
        hl2 = 1 - 1e-6  # f(D) >= 0
        half = 0
        if f(hl1) * f(hl2) < 0:
            while error > 1e-6:
                half = (hl2 + hl1) / 2  # 中间点
                val = f(half)
                if val < 0:
                    hl1 = hl1 + (hl2 - hl1) / 2
                elif val > 0:
                    hl2 = hl2 - (hl2 - hl1) / 2
                else:
                    return half
                error = np.abs(val)
            return half
        else:
            if self.alpha_L >= 0.24:
                return 1
            else:
                return 0

            # print("无法求得持液率稳态解, 使用其他方法进行流型判断")
            # return -1

    def judge_circle_flow(self) -> str:
        """判断环状流向段塞流得转换, 至于是不是段塞流还需要后续判定

        Returns:
            str: "AnnularFlow" 和 "non-AnnularFlow"
        """
        X2 = self.Lockhart_Martinelli2()
        Y = self.dimless_gravity()
        alpha_L = self.steady_alpha_L()
        if alpha_L == -1:
            if self.alpha_L >= 0.24:
                return "non-AnnularFlow"
            else:
                return "AnnularFlow"
        else:
            judge = Y * alpha_L**3 - (2 - 1.5 * alpha_L) / (1 - 1.5 * alpha_L) / X2
            if (judge > 0) or (self.alpha_L >= 0.24):
                return "non-AnnularFlow"
            else:
                return "AnnularFlow"

    #!==================================================================
    #!                       段塞流向弥散泡状流转换
    #!==================================================================
    def mixed_velocity(self) -> float:
        """计算混合速度

        Returns:
            float: 混合速度
        """
        return self.discount_V_g() + self.discount_V_l()

    def mixed_density(self) -> float:
        """计算混合密度

        Returns:
            float: 混合密度
        """
        return self.rho_L * self.alpha_L + self.rho_g * self.alpha_g

    def mixed_mu(self) -> float:
        """计算混合粘性

        Returns:
            float: 混合粘性
        """
        return self.mu_L * self.alpha_L + self.mu_g * self.alpha_g

    def mixed_re(self) -> float:
        """计算混合雷诺数

        Returns:
            float: 混合雷诺数
        """
        return self.mixed_density() * self.mixed_velocity() * self.D / self.mixed_mu()

    def mixed_friction_coefficient(self) -> float:
        """计算混合摩擦力系数

        Returns:
            float: 混合摩擦力系数
        """
        rhs = -2 * np.log(
            self.eps / 3.7065 / self.D
            - 5.0452
            / self.mixed_re()
            * np.log(
                1 / 2.8257 * (self.eps / self.D) ** 1.1098
                + 5.8506 / self.mixed_re() ** 0.8981
            )
        )
        return (1 / rhs) ** 2

    def cal_bubble_dmax(self) -> float:
        """计算段塞弥散为小气泡的最大气泡直径

        Returns:
            float: 最大气泡直径
        """
        dmax = (
            (4.15 * (self.discount_V_g() / self.mixed_velocity()) ** 0.5 + 0.725)
            * (self.sigma / self.rho_L) ** 0.6
            * (
                2
                * self.mixed_friction_coefficient()
                * self.mixed_velocity() ** 3
                / self.D
            )
            ** (-0.4)
        )
        return dmax

    def cal_dcd(self) -> float:
        """小于该直径的气泡不会变形，从而避免聚集与合并长大

        Returns:
            float: 该临界直径
        """
        return 2 * (0.4 * self.sigma / (self.rho_L - self.rho_g) / self.g) ** 0.5

    def cal_dcb(self) -> float:
        """由于浮力作用，大于该直径的气泡将向管子上部移动并形成分层从而向段塞流转换

        Returns:
            float: 该临界直径
        """
        return (
            3
            / 8
            * self.rho_L
            * self.mixed_friction_coefficient()
            * self.mixed_velocity() ** 2
            / (self.rho_L - self.rho_g)
            / self.g
            / np.cos(self.theta())
        )

    def judge_bubble_flow(self) -> str:
        """段塞流向弥散泡状流转换

        Returns:
            str: "BubbleFlow" 或者 "SlugFlow"
        """
        dmax = self.cal_bubble_dmax()
        dcd = self.cal_dcd()
        dcb = self.cal_dcb()
        if self.alpha_g <= 0.52:
            if (self.theta() <= 10 / 180 * np.pi) or (
                self.theta() >= 2 * np.pi - 10 / 180 * np.pi
            ):
                if (dmax < dcd) and (dmax < dcb):
                    return "BubbleFlow"
                else:
                    return "SlugFlow"
            else:
                return "SlugFlow"
        else:
            return "SlugFlow"

    # =====================================================
    #                   最终判断流型的函数
    # =====================================================
    def judge_pattern(self) -> str:
        """流型判断主程序

        Returns:
            str: "StratifiedFlow", "AnnularFlow", "BubbleFlow" 或者 "SlugFlow
        """
        progress1 = self.judge_layered_flow()
        if progress1 == "non-StratifiedFlow":
            progress2 = self.judge_circle_flow()
            if progress2 == "non-AnnularFlow":
                return self.judge_bubble_flow()
            else:
                return progress2
        else:
            return progress1
